Investigating the Impact of Irrigation on Malaria Vector Larval Habitats and Transmission Using a Hydrology‐Based Model

Abstract A combination of accelerated population growth and severe droughts has created pressure on food security and driven the development of irrigation schemes across sub‐Saharan Africa. Irrigation has been associated with increased malaria risk, but risk prediction remains difficult due to the heterogeneity of irrigation and the environment. While investigating transmission dynamics is helpful, malaria models cannot be applied directly in irrigated regions as they typically rely only on rainfall as a source of water to quantify larval habitats. By coupling a hydrologic model with an agent‐based malaria model for a sugarcane plantation site in Arjo, Ethiopia, we demonstrated how incorporating hydrologic processes to estimate larval habitats can affect malaria transmission. Using the coupled model, we then examined the impact of an existing irrigation scheme on malaria transmission dynamics. The inclusion of hydrologic processes increased the variability of larval habitat area by around two‐fold and resulted in reduction in malaria transmission by 60%. In addition, irrigation increased all habitat types in the dry season by up to 7.4 times. It converted temporary and semi‐permanent habitats to permanent habitats during the rainy season, which grew by about 24%. Consequently, malaria transmission was sustained all‐year round and intensified during the main transmission season, with the peak shifted forward by around 1 month. Lastly, we evaluated the spatiotemporal distribution of adult vectors under the effect of irrigation by resolving habitat heterogeneity. These findings could help larval source management by identifying transmission hotspots and prioritizing resources for malaria elimination planning.


Introduction
This file contains additional information on data collection, model development, calibration and simulation results.

Text S1. ParFlow-CLM Overview
ParFlow-CLM is a physical-based distributed model that couples an integrated hydrologic model, ParFlow, with a land surface model, Common Land Model (CLM).ParFlow simulates surface and subsurface flow in an integrated fashon by simultaneously solving Richards equation and shallow water equations (Kuffour et al., 2020) in three spatial dimensions.This is done by using Richard's equation to simulate variably saturated subsurface flow with a free surface overland flow boundary condition to accurately account for the surface-subsurface interactions (Kollet & Maxwell, 2006).The overland flow boundary condition can be activated and deactivated whenever the surface layer becomes saturated or unsaturated.Both saturation excess and infiltration excess mechanisms are allowed and streams occur naturally without the need to predefine river locations.In this study, the chosen boundary condition combines the diffusive wave approximation of the shallow water equations and Mannings equation.All fluxes simulated In ParFlow are driven by gradients in the hydraulic head.
CLM solves the land-water-energy balance, which includes evaporation, transpiration, snow processes, heat fluxes, and radiation partitioning (Maxwell & Miller, 2005).The two models can be connected across user-defined soil layers.Evapotranspiration and root uptake fluxes are calculated by CLM and incorporated in ParFlow as water fluxes into and out of the model through source or sink terms in Richards equation for subsurface flow.The fluxes are influenced by land surface and vegetation type which is characterized by a suite of parameters including soil color, canopy roughness, and green leaf area.
More details on ParFlow-CLM can be found in Kuffour et al. (2020).

Text S2. EMOD Malaria Model Overview
EMOD is an agent-based, mechanistic model which tracks interacting individuals in the form of humans and vectors and their environment to simulate malaria transmission at the population level (Bershteyn et al., 2018).This differs from a compartmental model which classifies humans and vectors into different homogenous population groups referred to as compartments and tracks the population based on a system of differential equations.
In EMOD each simulated individual goes through the various infection states of the susceptible-exposed-infected-recovered-susceptible epidemiological model (SEIRS).During each discrete 1-day time step of the simulation, the vector population evolves through maturation, mating, feeding and death within the model domain referred to as a node.
Successful feeding results in human infection and in turn, the malaria parasite is transmitted from the infectious human to the next susceptible vector.The node can be configured to represent any desired geographic scale ranging from the household level to the country level and the simulation can also include multiple nodes, with migration of both vector and human population enabled between nodes (Eckhoff, 2011;Eckhoff & Wenger, 2016).
The model is built upon six components namely, climate, larval habitat, vector transmission model, malaria infection and immune model, malaria symptoms and diagnostics and intervention (Bill & Melinda Gates Foundation, 2023).Climate can be configured for each node based on pre-set classification system tied to the native vegetation or user uploaded data.
In our case, we used climate data sourced from PERSIANN CCS-CDR (Sadeghi et al., 2021).In the default EMOD habitat model, there are temporary, semi-permanent and permanent natural habitats which are influenced by climate through rainfall and temperature and human -driven habitats which are a function of human population.The larval habitat area determines the number of eggs that develop into larvae as well as the number of emerging adult vectors.In the vector transmission model, the vector population is tracked throughout the mosquito life -cycle from egg to larvae and finally adult, which then progresses to the feeding cycle until death, comprising host-seeking, feeding, and egg-laying (Bill & Melinda Gates Foundation, 2023a, 2023b).The malaria infection and immune model tracks the parasite count for each infected individual, taking into account the immune response to infection.The model also tracks malaria symptoms such as fever and anemia and parasite count from diagnostics test to evaluate disease severity and mortality (Bill & Melinda Gates Foundation, 2023c).A clinical case is established when fever surpasses a certain threshold.Lastly, the intervention component allows the effect of antimalarial drugs, larviciding, long-lasting insecticidal nets (LLINs) and indoor residual spraying (IRS) to be modelled.For example, users can configure larviciding in the model based on a larval killing efficacy which is then applied on the larvae population through a decay function.

Text S3. Larval density estimation
In the study area, the surveyed larval habitats include drainage ditch, river edge/reservoir shoreline, swamp/marsh, rice puddle, animal footprint, tire track/road puddle, man-made pond, natural pond/rain pool, rock pool, water container, irrigation canal, and brick pit.The larval habitats were classified as temporary, semi-permanent, or permanent based on their natural characteristics.Since larval density can be significantly different in the dry and rainy season s (Hinne et al., 2021;Kweka et al., 2012) and the timing and duration of the survey periods were inconsistent, we sorted the measured larval densities from the 769 sample points (Figure 1) into the dry season (January to April; November to December) and the rainy season (May to October).We then calculated the average larval density for each season as shown in Figure S4.
In the surveyed area, the larval density for temporary habitats was higher in the rainy season than in the dry season during which the habitats are less stable.On the other hand, the larval densities for semi-permanent and permanent habitats were higher in the dry season.Most of the semi-permanent and permanent habitats were associated with river edges and swamps, whereby the larvae are prone to flushing in the rainy season.Finally, the larval density in Table 1 was calculated based on the average dry season and rainy season densities.

Text S4. Irrigation schedule design
Figure S6 shows the monthly irrigation schedule obtained from the Arjo-Didessa Sugar Factory, which is tailored to the sugarcane planting cycle.
To model irrigation in ParFlow-CLM, the irrigation interval and rate are required user inputs.A design report provided by the factory recommended 8-12 days for the design of the local irrigation system, so we set the irrigation interval as 10 days.
To determine the irrigation rate, we first calculated the irrigation depth, defined as the amount of water that needs to be applied when the soil water content is depleted to the wilting point.The irrigation depth () was calculated as where  is the field capacity,  is the permanent wilting point and ℎ  is the soil depth.
The study area is characterized by clay and clay loam with low permeability.Based on resources by the Northeast Region Certified Crop Advisor (Cornell University, 2010), the field capacity volumetric soil moisture content of clay was set as 50%, and the wilting point volumetric soil moisture content was set as 15%.A soil depth of 2 m was assumed.Using Equation (S1), an irrigation depth of 700 mm was obtained.
We configured the irrigation to be applied when 50% of the irrigation depth was depleted; hence, the actual irrigation depth to be applied over the 10-day irrigation interval was 350 mm.Adopting an intermittent irrigation strategy, we set the irrigation to be applied for 22 hours a day over 3 days within the 10-day cycle.Each irrigation period was assumed to be 3 days in order to match the local design irrigation rate of 5.1 mm/hr.Based on this assumption, the irrigation rate was calculated to be 5.3 mm/hr.

Text S5. Derivation of scaling factor and decay rate in the EMOD default habitat equations
In EMOD, the scaling factor (  ,   and   ) and decay rate (  and   ) are usually calibrated to match field observations of the total area of all larval habitats.As it is challenging to obtain this data for the study area, we used the simulated habitat area from Integrated EMOD to calibrate scaling factor and decay rate.As the permanent habitat area in EMOD is constant,   is set as the temporal mean of the fractional study area covered by permanent habitat (
For   and   , the default values in EMOD are 0.05 and 0.01.We conducted a sensitivity analysis of the intra-annual and inter-annual variability of fractional habitat area for   and   ranging from 0.01 to 0.2 as shown in Figure S7.For each iteration,   and   were set so that the mean temporary/semi-permanent habitat area matched ParFlow-CLM.Intra-annual variability was calculated in terms of the standard deviation of the 20-year average habitat area for each day of the year.Inter-annual variability was characterized by the standard deviation of the annual average habitat area for each year.Based on the sensitivity analysis results,   was set as 0.07 to match both the intra-annual variability and inter-annual variability of habitats from ParFlow-CLM as closely as possible.  was set as 0.01 by the same logic but it also has to be smaller than   since semi-permanent habitats should conceptually decay slower than temporary habitats.

Text S6. Model calibration
At each grid cell, ponding is assumed to occur if the soil saturation exceeds the threshold, .Therefore, the threshold was calibrated to ensure that the model will predict the occurrence of ponding at locations in line with the field-surveyed larval habitats for soil saturation above .The value for  was obtained based on a sensitivity analysis by altering the threshold and noting the corresponding change in the probability of detection ( ).The  determines if the model can predict an aquatic habitat successfully and can be calculated based on the ratio of the number of successful predictions or hits, , to the total number of samples, : Figure S8 shows the results of the sensitivity analysis.Generally, the  curve is higher for the simulation excluding dry season.This is because irrigation was only approximated by a simplified scheme in the dry season and may not reflect the localized irrigation dynamics.As the threshold was lowered,  increased because ponding occurred across a larger area in the model.The influence of topography on the ponding was weakened, and the soil type became the dominant factor.On the other hand, when the threshold was increased, less ponding was predicted, resulting in a lower  but the topographic variability was better represented.
Therefore, we selected a threshold of 0.75 for a reasonable POD of 0.66 (excluding dry season) without obscuring topographic variability.
In EMOD, we calibrated 15 key parameters identified from a preliminary sensitivity analysis, and Table S3 presents the calibrated values.Using the calibrated parameters, we compared the simulated prevalence rate against field data for Jan uary 2018 and October 2018 (Figure S9).The results are within the same order of magnitude.In addition, we compared the simulated monthly number of clinical cases with the recorded malaria cases from April 2018 to May 2020 (Figure S10).Apart from the two peaks missed in October 2018 and November 2019, the simulated malaria cases compare reasonably well with observation in terms of magnitude and pattern.As the clinical malaria cases were sourced from major hospitals within the study area, the two peaks in recorded cases could be anomalous due to an influx of patients from outside seeking treatment at the hospitals.Overall, the model shows good agreement with the field observation.

Text S7. Habitat Seasonality and Implications on Transmission
To evaluate the effect of the degree of seasonality in the larval habitat on malaria transmission, we conducted a sensitivity analysis using a synthetic sinusoidal time series for larval habitat fractional area with the same mean but different amplitudes: where  is the amplitude of fractional area,    is the fractional area at time ,  ̅ is the mean fractional area specific to the study derived from the hydrologic model.
The sensitivity analysis results can be found in Figure S14 and are summarized in Table S5.By reducing the seasonal amplitude from 0.2 to 0.1, the adult vector population remained relatively unchanged, but the vector infection and prevalence rates tripled.For the extreme case when  was reduced to 0, the vector infection and prevalence rates increased further by 4.29 times and 4.80 times, respectively.This finding agrees with the higher simulated malaria transmission in Default EMOD compared to Integrated EMOD (Figure S11e-g).It is possible that in the case where  was 0, the consistent adult vector population arising from the invariant habitat availability resulted in a stable parasite transmission throughout the year.As  increased, the disparity between the high and low vector abundance seasons increased.In the low vector abundance season, the transmission was minimal.In the high vector abundance season, transmission increased but would be limited by the human population.This resulted in an overall lower annual average vector infection rate and prevalence.Therefore, a nuanced approach considering the trend of the mean and the degree of the seasonality of larval habitats is required to predict malaria transmission accurately.Soil Type 250-meter -SoilGrids250m, TAXOUSDA (Hengl et al., 2017) Depth to Bedrock 250-meter -SoilGrids250m, BDRICM (Hengl et al., 2017) Near Surface Permeability (< 100 m) Regional Scale -GLobal HYdrogeology MaPS 2.0 (GLHYMPS, 2.0) (Gleeson et al., 2014) Table S2.Field data for ParFlow-CLM and EMOD validation.

Parameter Value
Antibody

Figure S1 .
Figure S1.Annual climate data from PERSIANN-CCS-CDR and ERA5 for the study area.(a) total precipitation, (b) average temperature and (c) average relative humidity.The red dashed line represents the linear trendline in each subplot.

Figure S2 .
Figure S2.Monthly climate data (averaged from 1994 to 2020) derived from PERSIANN-CCS-CDR and ERA5 climate data for the study area.(a) total precipitation, (b) average temperature, and (c) average relative humidity.

Figure S4 .
Figure S4.Percentage distribution of International Geosphere-Biosphere Programme type from land use survey.

Figure S5 .
Figure S5.Average larval densities for temporary, semi-permanent, and permanent habitats during rainy and dry seasons from field survey.

Figure S7 .
Figure S7.Adjustment of decay factors for temporary (  ) and semi-permanent habitats (  ) in default EMOD function based on sensitivity of (a) intra-annual variability and (b) interannual variability of habitats.

Figure S8 .
Figure S8.Sensitivity analysis of the probability of detection to saturation threshold.The probability of detection determines if the model can predict an aquatic habitat successfully and can be calculated based on the ratio of successful predictions to the total number of observations.The dotted vertical line corresponds to the selected threshold of 0.75, which results in a reasonable POD of 0.66 excluding dry season and a POD of 0.59 including dry season.

Figure S9 .
Figure S9.Comparison of simulated monthly average prevalence rate in Irrigation and measured prevalence diagnosed by Polymerase Chain Reaction (PCR).The whisker on the bar plot represents one standard error.

Figure S10 .
Figure S10.Comparison of simulated monthly confirmed cases in Irrigation and clinical data.The orange band indicates the 95% confidence interval.

Figure S11 .
Figure S11.Time series of daily climate data and comparison of simulated daily malaria transmission results between Default EMOD and Integrated EMOD.We only show the truncated time series from 2010 to 2020 instead of the full simulation period for brevity and readability.Climate data include (a) precipitation and (b) temperature.Malaria transmission results include (c) habitat larval capacity, (d) adult vector abundance, (e) adult vector infection rate, (f) entomological inoculation rate and (g) parasite prevalence rate.

Figure S12 .
Figure S12.The distribution of the top two-meter soil types in USDA soil taxonomy from SoilGrids250m TAXOUSDA dataset.Most soil types in this area are characterized as clay or clay loam with low permeability ranging from 0.0015 to 0.015 m/h.

Figure S13 .
Figure S13.Comparison of simulated semi-permanent habitats between Default EMOD (black line) and Integrated EMOD (blue line) in 2018.Earlier rising limb in Default EMOD: (1) no infiltration, new ponds created instantaneously by rainfall; (2) ponds formed sometime after rainfall when soil saturation exceeded the threshold.Delayed falling limb in Default EMOD: (3) habitat area continued to increase with rainfall; (4) pond drained/dried up and soil became unsaturated after a period without rainfall, new rainfall insufficient to create ponding as soil saturation remained below threshold.

Figure S14 .
Figure S14.Simulation results from sensitivity analysis of malaria transmission to different amplitudes of larval habitat seasonality, .Time series include (a) daily temperature, (b) synthetic sinusoidal larval habitat, (c) habitat larval capacity, (d) adult vector abundance, (e) adult vector infection rate and (f) parasite prevalence rate.

Table S4 .
Spatial average of adult vector from the dry season (November 2016 to April 2017) and the rainy season (May 2017 to October 2017).

Table S5 .
Average simulated adult vector abundance, adult vector infection rate, and parasite prevalence rate for different amplitudes of larval habitat seasonality, .